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ABSTRACT 

Context. We investigate the hydrodynamics of the core helium flash near its peak. Past research concerned with the dynamics of this 
event is inconclusive. However, the most recent multidimensional hydrodynamic studies suggest a quiescent behavior and seem to 
rule out an explosive scenario. 

Aims. Previous work indicated, that depending on initial conditions, employed turbulence models, grid resolution, and dimensionality 
of the simulation, the core helium flash leads either to the disruption of a low-mass star or to a quiescent quasi-hydrostatic evolution. 
We try to clarify this issue by simulating the evolution with advanced numerical methods and detailed microphysics. 
Methods. Assuming spherical or axial symmetry, we simulate the evolution of the helium core of a 1.25M Q star with a metallicity 
Z=0.02 during the core helium flash at its peak with a grid-based hydrodynamics code. 

Results. We find that the core helium flash neither rips the star apart, nor that it significantly alters its structure, as convection plays 
a crucial role in keeping the star in hydrostatic equilibrium. In addition, our simulations show the presence of overshooting, which 
implies new predictions concerning mixing of chemical species in red giants. 
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1. Introduction 

In stars of mass 0.7 M < M < 2.2 M Q the onset of helium burn- 
ing constitutes a major event - the core helium flash. The pre- 
flash stellar core contains a white dwarf-like degenerate struc- 
ture with a central density of about 10 6 g cirr 3 , and an off- 
center temperature maximum resulting from plasma- and photo- 
neutrino cooling. When helium burning commences in this de- 
generate core, the liberated nuclear energy cannot be used to ex- 
pand and cool the layers near the temperature maximum. Instead 
it causes further heating and a strong increase of the nuclear en- 
ergy release. Only when convection sets in, part of the excess 
energy can be transported away from the burning regions, in- 
hibiting thereby a thermonuclear explosion. At the end of the 
flash, the core has been expanded to densities of the order of 
10 4 g cm" 3 , with helium burning quiescently in the center, and 
the star has settled on the horizontal branch. While standard stel- 
lar evolution calculations have been very successful in reproduc- 
ing observations of stars on the main sequence and the red giant 
branch (RGB), we are forced to recognize several discrepancies 
concerning the post-flash phases. In particular, we recall the lack 
of understanding of the horizontal-branch morphology, of low- 
luminosity carbon stars, and of hydrogen-deficient stars. Since 
all these (and other) problems appear after the RGB phase, it is 
plausible to suspect that the helium flash may be treated incor- 
rectly in standard (hydrostatic) stellar evolution calculations. 

The conceptual problems associated with the helium core 
flash arise from the extremely short timescales involved in 
the event. While the pre-fiash evolution proceeds on a nuclear 
timescale of ~10 8 yrs, typical e-folding times for the energy 
release from helium burning can become as short as hours at 
the peak of the flash. Such short times are comparable to con- 



vective turnover times, i.e., the common assumptions used for 
the treatment of convection in stellar evolution codes (instanta- 
neous mixing, time-independence) are no longer valid. In addi- 
tion, the assumption of hydrostatic equilibrium no longer needs 
to be fulfilled. Early attempts to overcome these assumptions 
by allowing for one-dimensional hydrodynamic flow ( [Edwards] 
1 1 969} |Zimmermann| 1 970) | Villere| 1 976| | Wicketfj 1 977| > remained 
inconclusive. The results ranged from a confirmation of the stan- 
dard picture to a complete disruption of the star. 

|Cole & Deupree (T980, 1981 ) performed a two-dimensional 
hydrodynamic study of the core helium flash. However, their 
study was limited by the computational resources available at 
that time to a rather coarse computational grid (23 x 4 zones), 
a diffusive first-order difference scheme (weighted donor cell), 
and a short time evolution (10 5 s compared to the duration of 
the core helium flash of 10 s from the onset of convection). 
They observed, at the radius of the off-center temperature max- 
imum, a series of thermonuclear runaways where heat transport 
by convection and conduction was sufficiently efficient to limit 
the rise of temperature. Each runaway modified the convective 
flow pattern and led to some inward transport of heat across the 
off-center temperature inversion. During the simulation the time 
interval between runaways continuously shortened, and the max- 
imum temperature steadily increased until it eventually exceeded 

Deu pree & Cole| ([1983) and ( |Deupree1| 1 984a|b[ ) confirmed 
these findings using two-dimensional models with an improved 
angular resolution (6° instead of 20°), and three-dimensional 
simulations (with 8x8 angular zones in a 80° x80° cone, i.e., 10° 
angular resolution). Cole et aT] ( |1985| l performed stellar evo- 
lution calculations of the core helium flash using a model for 
convective overshooting based on these hydrodynamic Simula- 
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Fig. 1. Theoretical evolutionary track of a 1.25 M star with a 
metallicity Z=0.02 in the H-R diagram. The core helium flash 
begins at the tip of the red giant branch indicated by the arrow. 



tions. They found that the evolution of the core helium flash 
is unchanged except for about the last week prior to its peak. 
Furthermore, the possibility of mixing of core material into the 
hydrogen shell was suggested by numerical experiments where 
point source explosions were enforced (Deupree 1984b 1986; 
Deupree & Wallace 1987). These results raised the hope that 
some problems concerning abundance anomalies and mass loss 
could be solved by understanding the core helium flash. 

The results of the hydrodynamic simulations, though vary- 
ing in details, indicated a dynamic flash that could disrupt the 
star ( Deupree|1984a[) or at least lead to a significant loss of the 
envelope ( Cole & Deupree|1981|l. The simulatio ns were critized 
by Iben & Renzini (1984) and Fujimoto et al. (1990) because 
(i) the radial grid was too coarse, (li) the gravitational poten- 
tial was "frozen in" {i.e., time-independent), and (iii) because a 
"closed" outer boundary was used. The latter two assumptions 
tend to underestimate the expansion of the core, and hence tend 
to overestimate the violence of the flash. 

Since the work of Deupree the computational capabilities 
have grown tremendously and methods to simulate hydrody- 
namic flow have improved considerably. Thus, the limitations of 
the early studies concerning the grid resolution and the numeri- 
cal treatment, which were the main points of critique, meanwhile 
can be reduced considerably. At the same time, we still have no 
coherent picture up to what extent and under what circumstances 
(stellar mass and composition) hydrodynamic core helium flash 
evolution could differ from canonical stellar evolution calcula- 
tions. It therefore appears necessary to have a new and fresh look 



into the dynamics of the core helium flash. Incidentally, Deupree 
( |1996| ) himself re-examined the problem already more than a 
decade ago concluding that the flash does not lead to any hydro- 
dynamic event. Quiescent behavior of the core helium flash is 
also favored by recent three-dimensional simulations (Dearborn 
et al. 2006[ Lattanzio et al. 2006 1 where the energy transport due 
to convection, heat conduction, and radiation seems to be always 
able to transport most of the energy generated during the flash 
quiescently from the stellar interior to the outer stellar layers, 
implying no hydrodynamic event, and hence a quasi-hydrostatic 
evolution. 



Fig. 2. Temperature distribution as a function of radius. The 
dashed line gives the distribution obtained from stellar evolu- 
tionary calculations with the "Garstec" code, while the solid line 
shows the mapped and stabilized distribution used as initial con- 
dition in the hydrodynamic simulations. CVZ marks the convec- 
tion zone. 



In the following we present a completely independent inves- 
tigation of the core helium flash by means of one-dimensional 
and two-dimensional hydrodynamic simulations using state-of- 
the-art numerical techniques, a detailed equation of state, and a 
time-dependent gravitational potential. The hydrodynamic cal- 
culations cover about 8 hrs of the evolution near the peak of the 
core helium flash. In passing we note that the present investiga- 
tion was instigated by a similar, meanwhile technically obsolete 
study which was performed by Kurt Achatz ( Achatz 1995 ) in the 
context of his diploma thesis. The results of this latter study have 
unfortunately never been published. 

The paper is organized as follows. In Sect. 2 we discuss 
briefly the stellar input model for the simulations along with 
some results from hydrostatic core helium flash calculations. 
In Sect. 3 the hydrodynamics code and the numerical methods 
are introduced, while the results of our one and two-dimensional 
hydrodynamic runs are presented in Sect. 4 and 5, respectively. 
Finally, the conclusions are given in Sect. 6. 

2. Initial stellar models and hydrostatic calculations 

Table[T] summarizes some properties of our initial model, which 
was obtained from stellar evolutionary calculations with the 
"Garstec" code ( |Weiss & Schlattl||200"0"l [2007] ). It corresponds 
to a star with a mass of 1 .25 M Q and a metallicity Z = 0.02 at 
the peak of the core helium flash {Lfj e ~ 1O 9 L ) evolved with a 
hydrostatic stellar evolution code. During this violent episode, 
the star is located at the tip of the red giant branch in the H-R 
diagram (Fig.[T]), hence being a red giant consisting of a small 
central helium core with a radius r ~ 1.9 10 9 cm, surrounded by 
a hydrogen burning shell and a huge convective envelope with 
a radius i~ 10 13 cm. Figure[2] shows the temperature distribu- 
tion inside the helium core, which is characterized by an off- 
center temperature maximum T max , from where the temperature 
steeply drops towards smaller radii and follows a super-adiabatic 
gradient towards larger radii (convection zone). The radius r max 
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Table 1. Some properties of the initial model: total mass M, stellar population, metal content Z, mass Mhc and radius Rn e of the 
helium core (X( 4 He) > 0.98), nuclear energy production in the helium core Ln e , maximum temperature of the star T max , and radius 
r m ax and density p ma x at the temperature maximum. 
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Fig. 3. Left panel: Pressure (in 10 22 dyn crrT 2 ) and density (in 10 5 g cm -3 ) distribution of the mapped and stabilized initial model. 
The pressure and density profile of the original stellar evolution cannot be distinguished from the profiles of the mapped model on 
this scale. Right panel: Chemical composition of the initial model showing a dominant fraction of helium and an apparent peak in 
12 C at the position of the temperature maximum resulting from a non-instantaneous treatment of convective mixing. 



of the temperature maximum coincides with the bottom of the 
convection zone. The almost discontinuous temperature stratifi- 
cation near T max (temperature inversion), where the temperature 
rises from 7 10 7 K to 1.7 10 8 K, results from an interplay between 
neutrino cooling and heating by nuclear burning. Figure[3]shows 
the density and pressure stratification of the model. One recog- 
nizes that the temperature inversion is correlated with a drop in 
density. A detailed view reveals that the steep increase of tem- 
perature corresponds to a decrease of the density by 11%, an 
increase of the ion pressure by 70%, and a drop of the electron 
pressure by 9%, respectively. Even at the peak of the core he- 
lium flash, the helium core is still strongly degenerate: compared 
to the electron pressure the ion pressure is lower by a factor of 
6, while the radiation pressure is smaller by almost 3 orders of 
magnitude. 

The stellar model contains the chemical species 'H, 3 He, 
4 He, 12 C, 13 C, 14 N, 15 N, 16 , I7 0, 24 Mg, and 28 Si. However, 
since we are here not interested in the detailed chemical evolu- 
tion of the star, it is not necessary to consider all of these species 
in our hydrodynamic simulations, as the triple-a reaction dom- 
inates the energy production rate during the core helium flash. 
For our hydrodynamic simulations we thus adopt only the abun- 
dances of 4 He, 12 C, and ls O. The remaining composition is as- 
sumed to be adequately represented by a gas with a mean molec- 
ular weight equal to that of 2() Ne (Fig.[3jl. 

The stellar evolutionary model is one-dimensional, hydro- 
static, and was computed on a Lagrangian grid of 2294 zones. 
Since only the helium core of the model (without its very cen- 
tral part; see Sect.|3.6|> is of interest to us, we consider only the 



initial data for 2 10 8 cm < r < 1.2 10 9 cm, and interpolate all 
relevant quantities {e.g., density, temperature, composition) onto 
our Eulerian, lower resolution computational grid using polyno- 
mial interpolation ( [Press et aL 1992| l. Due to the interpolation er- 
rors and subtle differences in the input physics, the interpolated 
model is no longer in perfect hydrostatic equilibrium. In order 
to perfectly balance also the gravitational and pressure forces 
in the interpolated model, we use an iterative procedure in the 
first hydrodynamic timestep to minimize the numerical fluxes 
across zone boundaries. The whole process results in a small 
temperature decrease with respect to the temperature profile of 
the original model (Fig.|2]i. The differences do not exceed a few 
percent depending on the radial resolution of the Eulerian grid. 
The resulting changes in the density and pressure profiles are 
negligible due to the strong electron degeneracy of the gas. The 
main cause for the slight de-stabilization of the mapped initial 
stellar model is the use of different equations of state in both 
codes. The hydrodynamic code employs the equation of state by 



Timmes & Swesty (2000i, whereas the "Garstec" code relies on 
the OPAL equation of state by Rogers et al.| ( |1996[ ). At a given 
density, temperature, and composition in the helium core during 
the flash, these equations of state give pressure values which dif- 
fer typically by 1 % the difference being most apparent in regions 
where the matter is highly degenerate. 

Given that the maximum temperature in the helium core is 
T ~ 1 10 H K, the stellar model reaches the peak in nuclear en- 
ergy production rate during the core helium flash in less than 
10 4 yrs. The rate at which the nuclear energy production rises is 
highly non-linear. From the onset of the core helium flash at a 
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3.2. Neutrino emission 
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Fig. 4. Temporal evolution of the helium luminosity L# e (solid) 
versus the hydrogen luminosity L# (dash-dotted) of model M 
during the core helium flash. 



helium luminosity Ln e ~ lO 1 /^, it takes almost 30000 yrs to 
reach Ln e ~ 10 4 L G , whereas it requires only 40 yrs to reach 
Lhc ~ 10 10 . The first core helium flash is followed by four 
subsequent mini flashes (Fig. [4]) identified as thermal pulses by 
|Thomas| (l967) until the degeneracy in the helium core is lifted 
completely and the star settles down on the horizontal branch 
quiescently burning helium in its core. 

Since the computed model is a Pop I metal rich star, it does 
not experience any hydrogen entrainment during the core helium 
flash ( [Fujimoto et : al.|1990| [SchTattl et al.|200l] >. 



3. Input physics and numerics 

3.1. Thermal transport 

The energy flux density due to thermal transport is given by 

/cond = -^condVr, (1) 

where K con< i is the total conductivity ( erg K -1 cirT 1 s~') and VT 
the temperature gradient. 

In the helium core, which is partially degenerate, thermal 
transport due to both radiative diffusion and electron conduction 
is important, while heat transport by ions is negligible, i.e., 



Ky — 



(2) 



(3) 



^cond — Ky + K e . 

The radiative conductivity is given by 

Aac T 3 

3 KyP ' 

where a:, a, and c are the Rosseland mean of the opacity, the radi- 
ation constant, and the speed of light, respectively. For the opac- 
ity, we use a fit formula due to|Iben| ([l975) which is based on 
the work by |Cox & Stewart] ( |1970b|a[ ). It takes into account the 
radiative opacity due to Thomson scattering, free-free (Krames 
opacity), bound-bound, and bound-free transitions. 

For the thermal transport by electron conduction we consider 
contributions due to electron-ion, and electron-electron colli- 



sions which are treated according to Yakovlev & Urpin ( 1980), 
and |Poteldiin"eraT1 ( fT997] l. 



The evolutionary time covered by our hydrodynamic simulations 
is too short for neutrino cooling to be of importance. The neu- 
trino losses computed from the analytic fits of |Itoh et al.| (l996) 
give a cooling rate e < 10 2 erg g -1 s _1 , or a corresponding de- 
crease of the maximum temperature by \tST\ < 10~' K over the 
longest simulations we performed. Hence, cooling by neutrinos 
was neglected. 



3.3. Equation of state 

The equation of state employed in our hydrodynamic code 
includes contributions due to radiation, ions, electrons, and 
positrons. Thus, the total pressure is given by 



5479.0 where 



P — Py P ion "t" P e *p > 



(4) 



(5) 



is the radiation pressure of a black body of temperature T (a is 
the universal radiation constant), and 



(6) 



is the pressure of a non-relativistic Boltzmann gas of density p 
consisting of a set of ions of abundance Y, = Xj/Aj (X, and A, 
are the mass fraction and the atomic mass number of species i, 
respectively). P e + P p is the pressure of an arbitrarily degenerate 
and relativistic electron-positron gas based on table interpolation 
of the Helmholtz free energy ( Timmes & Swestyj2 000). 



3.4. Nuclear burning 

The energy generation rate by nuclear burning is given by 

Ato,c 2 • 



z 



(7) 



where 

Am; = Mi - Aiin u . (8) 

is the mass excess of a nucleus of mass M,, and m u is the atomic 
mass unit. 

Abundance changes are described by a nuclear reaction net- 
work consisting of the four c-nuclei 4 He, 12 C, ie O, and 20 Ne, 
coupled by seven reactions (including the triple-o- reaction). We 
used the reaction rate library of Thielemann (private communi- 
cation), which gives the product of the Avogadro number and 
the velocity averaged cross section (crv) in terms of the fit for- 
mula 



N A {<rv) = J] exp[c u + c % T~ X + c 3 ,T- 1/3 + c 4 ,T 
+ c 5 ,T + c 6/ r 5/3 +c 7/ lnr] 



1/3 



1=1 



(9) 



with rate dependent coefficients c,y (1 < i < 7). Up to three sets 
of coefficients (i.e., 1 < «/ < 3) are used. The total reaction rate 
due to all one body, two body, and three body interactions has 
the form ( |Muller|1998| l: 



Yi-^cJ.Yj + ^iU^)pN A (crv}j tk YjY k 



Z 



c i (j,k,P,p 2 Nl(crv)j tk jYjY k Y l 



(10) 
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where the weight factors c,- inhibit multiple counts in the sums 
over the nuclei j,k,l. The following nuclear reactions were con- 
sidered: 



3.6. Code 
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Mathematically this results in a nuclear reaction network con- 
sisting of seven non-linear first order differential equations of 
the form given by Eq. ( p~0] > and a temperature equation 



dT 



dT 

7 de' 



(11) 



where s is the specific internal energy. 

The effects of electron screening were included according to 
Dewitt et al. ( 1973 1 for the triple-o- reaction rate, and in the weak 
screening regime only. 

3.5. Evolutionary equations 

The hydrodynamic and thermonuclear evolution of the core he- 
lium flash was computed by solving the governing set of fluid 
dynamic equations in spherical coordinates on an Eulerian grid. 
Using vector notation these equations have the form, 



3U 

~dt 



+ VF = S 



(12) 



with the state vector U 



U = 



P 
p\ 
pe 
pYi) 



(13) 



the flux vector F 



F = 



and the source vector S 



S = 



pv 
pvv 

(pe + p)y + f cond 
pYtV 





-pVO 
-pv ■ V0> + pe m 

PY 



(14) 



(15) 



with i — 1, . . . , A^nuc where N nuc is the number of nuclear species 
considered in the nuclear reaction network, and p, p, v and <1> are 
the density, pressure, velocity and gravitational potential, respec- 
tively. The term f com t describes energy transport by thermal con- 
duction (see Sect. 3.1), and e nuc and the F, are the nuclear energy 
generation rate and the change of the mass fraction of species 
i due to nuclear reactions, respectively (see Sect. 3.4). The total 
energy density pe - pe + pvv/2 with e being the specific total 
energy. 



The numerical simulations were performed with a modified ver- 
sion of the hydrodynamic code Herakles ( [Kifonidis et al.|2003~ 
2006), which is a descendant of the code Prometheus devel- 
oped by Bruce Fry xell and Ewald Miiller (Mii ller et al.||1991[ 
Fryxell et al.|1991 1. The hydrodynamic equations are integrated 



to second order accuracy in space and time using the dimen- 
sional splitting approach of [S trang ( 1968 ), the PPM reconstruc- 
tion scheme (Colella & Woodwar d|1984|l, and a Riemann solver 
for real gases according to |Colella & Glaz| ([r984). The evolution 
of the chemical species is described by a set of additional conti- 
nuity equations (Plewa & Miiller 1999). Source terms in the evo- 
lutionary equations due to self-gravity and nuclear burning are 
treated by means of operator splitting. Every source term is com- 
puted separately, and its effect is accounted for at the end of the 
integration step. The viscosity tensor is not taken into account 
explicitly, since the solution of the Euler equations with the 
PPM scheme corresponds to the use of a sub-grid scale model 
that reproduces the solution of the Navier-Stokes equations rea- 
sonably well (Meakin & Arnett [2007 1. Thermal transport is 



treated in a time-explicit fashion when integrating the evolution- 
ary equations. Self-gravity is implemented according to Miiller 
|& Steimnetz| (fl995), while the gravitational potential is approx- 
imated by a one-dimensional Newtonian potential which is ob- 
tained from the spherically averaged mass distribution. The nu- 
clear network is solved with the semi-implicit Bader-Deufelhard 
method which utilizes the Richardson extrapolation approach 



Press 



and sub-stepping techniques (Bader & Deuflhard 1983 
|et al.| 199*2] ) allowing for very large effective time steps. 

The code is vectorized and allows for an adjustment of the 
vector length to the memory architecture. Therefore, an optimal 
performance on both vector and super-scalar, cache-based ma- 
chines can be achieved. 

A program cycle consists of two hydrodynamic timesteps 
and proceeds as follows: 

1 . The hydrodynamic equations are integrated in r-direction (r- 
sweep) including the effects of heat conduction. The time 
averaged gravitational forces are computed, and the momen- 
tum and the total energy are updated to account for the grav- 
itational source terms. Subsequently, the equation of state is 
called to update the thermodynamic state due to the change 
of the total energy. 

2. Step (1) are repeated in ^-direction (0-sweep). 

3. The nuclear network is solved in all zones with significant 
nuclear burning (T > 10 8 K). Subsequently, the equation of 
state is called to update the pressure and the temperature. 

4. In the subsequent timestep the order of Step (1) and (2) is re- 
versed to guarantee second-order accuracy of the time inte- 
gration, and Step (3) is repeated with the updated quantities. 

5. The size of the timestep for the next cycle is determined. 

When using spherical coordinates, the CFL stability condi- 
tion on the timestep is most restrictive near the origin of the 
grid. However, inside a region beneath the off-center tempera- 
ture maximum there are no significant non-radial motions to be 
expected during the evolution of the core helium flash except 
in the immediate vicinity of the temperature inversion, where 
convective overshooting may occur. Hence, cutting out the very 
center of the computational grid does not lead to any numerical 
bias, but saves considerable amounts of computational time. In 
the radial direction we used a closed (i.e., reflective) outer and 
inner quasi-hydrostatic boundary obtained by means of polyno- 
mial extrapolation, which significantly suppresses any artificial 
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4. Results of 1D simulations 
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Time (10 3 s) 

Fig. 5. Evolution of the temperature maximum T max in the one- 
dimensional models JE2 (solid), JE3 (dashed), and JE4 (dash- 
dotted), respectively. 
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Fig. 6. Temperature stratification across the helium core in 
model JE4 during the runaway at t\ = 12270 s (dotted), ti = 
12352 s (dashed), f 3 = 12392 s (dash-dotted), and t 4 = 12762 s 
(dash-dot-dotted), respectively. The solid line corresponds to the 
initial model (to), and the arrow indicates the direction of the 
flame propagation. 



Table 2. Some properties of the ID simulations: number of ra- 
dial grid points (N r ), radial resolution (Ar in 10 8 cm), time up to 
the thermonuclear runaway, t trn , and maximum evolution time 
tmax (both in s). 



run 


N r 


Ar 




tmax 


JE2 


180 


5.55 


40700 


42500 


JE3 


270 


3.77 


14600 


16250 


JE4 


360 


2.77 


12300 


15600 



We have performed several one-dimensional simulations us- 
ing model M, which differ only by their grid resolution (see 
Table|2|l to see whether without allowing for convective flow a 
thermonuclear runaway can be avoided. 

Figure[5]demonstrates that heat conduction and adiabatic ex- 
pansion alone fail to stabilize the model, i.e., one-dimensional 
hydrodynamic simulations result in a thermonuclear runaway. 
Initially, the maximum temperature increases only slowly, but 
it starts to rise rapidly after a time t trn (Tab.|2]i up to a value 
T ~ 10 9 K. For instance, from the temperature evolution of 
model JE4 one can determine that a local hot spot with a tem- 
perature of 2.3 10 8 K will runaway after about 80 s (Fig. [6]). The 
time at which the runaway is triggered depends on the grid reso- 
lution, being longer in models with lower resolution (Fig.BJ. 

In every case, a thermonuclear flame with T ~ 10 9 K ulti- 
mately forms and propagates outwards with a subsonic veloc- 
ity depending on the grid resolution. Since our two-dimensional 
(more realistic) simulations do not show such a behavior, we will 
refrain from further discussing details of the one-dimensional 
simulations. 



5. Results of 2D simulations 

In Table|3]we summarize some characteristic parameters of our 
two-dimensional simulations that are based on model M. 

We will first discuss one specific simulation DV4 in some 
detail, which serves as a standard to which we will compare the 
results of other runs. Thereafter, we will discuss some general 
properties of all 2D simulations. Every simulation covered ap- 
proximately 30000 s (~ 8 hrs) of the evolution near the peak 
of the core helium flash. They were performed on an equidis- 
tant spherical grid encompassing 95% of the helium core's mass 
(X( 4 He)>0.98) except for a central region with a radius of 
r = 2 10 8 cm which was excised in order to allow for larger 
timesteps. As this radius is sufficiently smaller than the radius of 
the temperature inversion (r ~ 5 10 8 cm), its presence does not 
influence the convection zone. 



velocity fluctuations resulting from an imbalance of gravitational 
and pressure forces in the boundary zones. For two-dimensional 
runs, the boundary conditions in the angular direction are reflec- 
tive as well. 

After interpolation and stabilization, the initial model in the 
two-dimensional simulations had to be perturbed explicitly to 
trigger convection, because an initially exactly spherically sym- 
metric model remains that way for ever when evolved in spheri- 
cal coordinates with our code. We imposed a random flow field 
with a maximum (absolute) velocity of 10 cm s _1 , and random 
density perturbations with Ap/p < 10~ 2 . 



5.1. Simulation DV4 

After the start of the simulation the initial velocity perturbations 
begin to grow in a narrow layer just outside the temperature max- 
imum (r ~ 5 10 8 cm), i.e., in the region heated by nuclear burn- 
ing. Later on at t ~ 800 s, several hot bubbles appear, which 
rise upward with maximum velocities ~ 4 10 6 cm s _1 (Fig. [8}. 
They are typically about 0.2% hotter than the angular averaged 
temperature at a given radius. The 4 He mass fraction of all hot 
bubbles is about 0.4% less than the corresponding angular aver- 
aged value since helium has been depleted in the bubbles by the 
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Table 3. Some properties of the 2D simulations: number of grid points in radial (N r ) and angular (Ng), radial (Ar in 10 cm) and 
angular grid resolution (A8), characteristic length scale l c of the flow (in 10 8 cm), characteristic velocity v c of the flow (in 10 6 cm s _1 ), 
Reynolds number R n associated with the numerical viscosity of our code ( |Porter & Woodw ard 1994), damping time-scale due to 
the numerical viscosity t„, typical convective turnover time t a , and maximum evolution time t max (in s), respectively. 
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Fig. 7. Left panel: Temporal evolution of the horizontally averaged temperature maximum (T) max (solid), and of the global temper- 
ature maximum T max (dotted) in model DV4. The dashed line corresponds to the temporal evolution of the maximum temperature 
in the stellar evolutionary calculations of the model M. Right panel: The r.m.s convection velocity v cnv in simulation DV4 averaged 
over 6000 s (solid) versus the convection velocity predicted by the mixing length theory v,„/, (dashed). 



triple a reaction. Consequently, 12 C and 16 (produced in helium 
burning) are enhanced by ~ 0.7% in the bubbles. 

During the first 700 s of the evolution, the off-center maxi- 
mum mean temperature (T) max rises with a rate of ~ 1000 K s" 1 
until it reaches a value ~ 1.67 10 8 K. At this moment, from the 
region around the (T) max , the bubbles emerge and cause its de- 
crease by ~ 2.6 10 6 K in just 570 s corresponding to a tempera- 
ture drop rate of 4540 K s _1 (Fig. [7]). This phase marks the on- 
set of convection where a fraction of the thermonuclear energy 
released via helium burning starts to be efficiently transported 
away from the burning regions by matter flow, thereby inhibit- 
ing a thermonuclear runaway. 

Once the bubbles form, they rise upwards and start to inter- 
act and merge, i.e., the convective layer begins to grow in radius. 
About ~ 1300 s after the start of the simulation, the whole con- 
vection zone is covered by an almost stationary flow pattern with 
an almost constant total kinetic energy of the order of 10 45 erg. 
At this time vortices dominate the flow pattern. They extend 
across the whole convective region (~ 2AH p ), and are of ap- 
proximately similar angular size, one vortex covering about 40 
degrees (diameter ~ 5 10 8 cm). Usually we find about four such 
vortices with two dominant up-flows of hot gas at 8 ~ 60°, and 
9 ~ 120°, respectively (see, e.g., Fig. [8]). These large vortices are 
rather stable surviving until the end of our simulations. Typical 
convective flow velocities are v cm , ~ 1.5 10 6 cm s _1 , and thus 
well below the local sound speed (cs ~ 1.7 10 8 cm s~'), i.e., a 



vortex requires about 600 s for one rotation. The persistence of 
vortices is not typical for turbulent convection. 

The dominance of large scale structures might be a conse- 
quence of the usage of a Riemann solver based compressible 
code. The Mach number M of the convective flow is ~ 0.01. 
Is PPM suited for this kind of subsonic flow? This question, 
which is beyond the scope of the present study, needs to be in- 
vestigated, as it is know that the artificial viscosity of standard 
Riemann solver methods exhibit incorrect scaling with the flow 
Mach number as M — > 0. (Turkel 1999) i.e., the inherent artifi- 
cial viscosity of PPM may be too high for adequately simulating 
flows at low Mach numbers {e.g., M ~ 0.01). 

Energy transport by convection within the vortices is concen- 
trated into a few narrow upward drafts, compensated partially, 
but only to a small extent, by down-flows. The vortices transport 
energy mostly along their outer edges. Matter in their centers 
does not interact with regions of dominant nuclear energy pro- 
duction at all. 

The horizontally averaged value of the maximum temper- 
ature, barring some additional temperature fluctuations due to 
convection, is slightly rising after the onset of convection during 
the whole subsequent evolution with a rate of around 40 Ks" 1 
(see Fig.[7](. This rate seems to be about 60% smaller than the 
rate seen in the stellar evolutionary calculations (~ 100 K s _1 ), 
which could be either a result of the initially lower value of the 
temperature maximum after the stabilization phase at the begin- 
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Fig. 8. Snapshots of the onset of convection at 1020 s (upper panels), and of the evolved convection at 29000 s (lower panels) in 
model DV4, showing the temperature contrast AT = 100 (T - (T) g )/(T) g (left panels), the velocity field (middle panels), and the 
12 C contrast A 12 C = 100 ( 12 C - ( l2 C)g)/( 12 C)g (right panels), respectively. denotes a horizontal average at a given radius. 



ning of the simulation (see Sec. 2) or more dynamic convective 
motion, since the mean convective velocities v„„. exceed the ve- 
locities predicted by mixing length theory, v m /,, on average by a 
factor of four (Fig.|7jl. 

Convection distributes the energy in such a way that the tem- 
perature gradient V never significantly exceeds V fl </ in model M. 
Although, the value of V established at the beginning of the sim- 
ulation deviates slightly after some time from the gradient at 
later times, it remains close to the adiabatic temperature gradi- 
ent V fl£ / (the relative difference is less than 1 %). In this respect 



there is thus no indication of any significant deviation from the 
situation obtained in stellar evolutionary calculations. 

The apparent spike in the initial 12 C distribution at the loca- 
tion of the temperature maximum (Fig.[3]l is a result of a non- 
instantaneous treatment of the convective mixing in stellar evo- 
lutionary calculations. It turns out that a non-instantaneous treat- 
ment of mixing is not required during the core helium flash since 
simulation DV4 indicates that the spike gets smeared out imme- 
diately after convection is triggered. This implies that the as- 
sumption of instantaneous mixing is a good approximation lo- 
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Fig. 9. Snapshots of various energy fluxes and source terms in model DV4 (time averaged over 6000 s from t = 18000 s to t — 
24000 s): (a) convective flux Fc (solid), and the energy flux due to the thermal transport Fp (dash-dotted); (b) kinetic flux F% 
(solid), acoustic flux Fp (dash-dot-dotted), and sum of the kinetic and convective flux Fc + Fk (dashed); (c) source terms due to 
work done by buoyancy forces Pa, and (d) due to volume changes Pp. The vertical lines enclose the nuclear burning zone (T > 10 s 
K). 



cally, despite the strong temperature dependence of the energy 
production rate. 

5.1 .1 . Energy fluxes 

Fig. [9] displays the individual contributions of various energy 
fluxes, time-averaged over many convective turnover times, 
i.e., only the average effect of convection should be apparent. 
The derivation of these quantities is explained in Appendix A. 
All energy fluxes, F, describe the amount of energy which is 
transported per unit of time across a sphere of given radius. 

Most of the nuclear energy production in the convection zone 
takes place in a relatively narrow shell around the location of 
the temperature maximum. This energy is transported away by 
both convection and thermal transport due to heat conduction 
and radiation. The convective (or enthalpy) flux, Fc, varies from 



-0.2 10 42 erg s" 1 up to 1.6 10 42 erg s" 1 . The kinetic flux, F K , 
reaches a value of at most 1 10 42 erg s _1 , and is mostly positive 
in the convection zone, i.e., the motion has a predominantly up- 
ward direction. This implies that the fast narrow upward directed 
streams are dominating over the slower and broader downward 
flows. The ratio of the extreme values of Fc and Fk is nearly 
2:1, i.e., nuclear energy is mainly stored in the internal energy of 
rising hot gas. 

Convective and kinetic energy flux together transport more 
than 90% of the generated nuclear energy upward through the 
convection zone, the value is dropping to zero towards its border. 
Part of the heat released in the nuclear processes is in fact trans- 
ported downwards towards the inner edge of the temperature in- 
version. Almost none of the nuclear energy reaches the surface 
of the helium core, neither by convection nor by conduction, 
i.e., all the energy released is deposited within the core causing 
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Fig. 10. Angular averaged 12 C distribution (dashed) and temperature stratification (thick) at the inner (left panel) and outer edge 
(right panel) of the convection zone in model DV4 at t — 30000 s. The vertical dotted lines mark the initial boundaries of the 
convection zone at t — s. 



its expansion. Energy transport due to heat conduction and radi- 
ation is everywhere negligible compared to the other contribu- 
tions. The viscous flux, Fy, is small as well, and losses due to 
friction, Py, only influence the dynamics significantly near the 
borders of the convection zone (Achatz 1995]). 

For completeness we also consider the flux and source terms 
of the kinetic energy (see Appendix A), which allow for a further 
insight into the operation of convection. The radial profile of the 
source term P A , corresponding to the work done by buoyancy 
forces, shows that the vertical convective flows are accelerated 
due to their density fluctuations in the entire region of dominant 
nuclear burning (burning zone) above T max . Corresponding pres- 
sure fluctuations (causing expansion due to a pressure excess, 
respectively compression due to a pressure deficit) powered by 
the volume work P P show that the gas within the burning region 
expands, which effectively again implies that an acceleration oc- 
curs. Due to the importance of Pp in the convection zone, the 
acoustic flux F P , which transports pressure fluctuations, reaches 
a value comparable to that of the kinetic flux Fg , its value being 
negligible elsewhere. 



5.1.2. Turbulent entrainment, temperature inversion and the 
growth of the convection zone 

Turbulent entrainment (commonly referred to as overshoot- 
ing) is a hydrodynamic process allowing for mixing and heat- 
ing in regions which are convectively stable according to 
the Schwarzschild or Ledoux criterium. Turbulent entrainment, 
i.e., penetration beyond the formal convective boundaries, takes 
place at both edges of the convection zone, and is driven by 
down-flows and up-flows. We study the entrainment by moni- 
toring the temperature changes and the 12 C concentration at the 
(formal) edges of the convection zone. I2 C is the most suitable 
element for investigating the extent of convective mixing, be- 
cause at the beginning of the simulations, it is mostly absent 
outside the convection zone, and therefore can be enhanced there 
only due to overshooting. 



At t — 30000 s, i.e., near the end of simulation DV4, the tern 



perature inversion is located at r — 4.65 10 8 cm (Fig. 1 10}. Thus, 
it is about 70 km closer to the center of the star than it was at the 
beginning of the simulation (4.72 10 8 cm). Its shape remains al- 
most unchanged and discontinuous during the whole evolution, 
and its propagation speed can be estimated from the heating rate 
6T/6t ~ 2760 K s _I and the local gradient ST/Sr ~ 12 K cnr 1 
at the steepest point of the inversion: 



-{ST I St) I {ST/Sr) ~ -2.3 m s" 



(16) 



This speed is significantly higher than the propagation speed due 
to the heat conduction alone. Note that the energy flux carried 
by the heat conduction is seven orders of magnitude smaller than 
the energy flux carried by the convection. Assuming that the con- 
vective energy flux at the position of the temperature inversion 
{F c ~ 0.2 10 42 erg s ) is used up completely to heat the lay- 
ers beneath the temperature inversion, a typical heating rate of 
t = E/Ci„ v ~ 1250 K s" 1 can be derived, which is a bit smaller 
than the value inferred from the simulation, but still in good 
agreement. C; nv is the heat capacity of the layers including the 
temperature inversion (C; nv ~ 1.6 10 38 ergKr 1 ). This implies 
that turbulent entrainment leads to a strong heating of the inner 
neutrino cooled center of the star that occurs on timescales which 
are relatively short compared to stellar evolutionary timescales. 
Such a heating was studied already by Deupree & Cole ( 1983) 



and |Cole et al.| ( fl985] l who obtained qualitatively similar results. 
Note, that in the one-dimensional stellar evolution calculations 
the temperature maximum moves outwards with time. 

Assuming that the estimated propagation speed of the tem- 
perature inversion remains constant, it would reach the center of 
the helium core and lift the electron degeneracy there in just 24 
days. This scenario would rule out the occurrence of mini-flashes 
subsequent to the main core helium flash, which are observed in 
stellar evolutionary calculations (Fig.|4|. Moreover, as in stars 
with higher mass and helium abundance the flash occurs closer 
to the center (Swei gart & Gross|19 78), in these stars the center 
can be reached even faster. 

We have also found an influence of the turbulent entrain- 
ment on the outer boundary of the convection zone. In the ini- 
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Fig. 11. Left panel: Evolution of the total energy production rate in solar luminosity L for models DV2 (dotted), DV3 (dashed), 
and DV4 (dash-dotted), respectively. Right panel: Mean temperature distribution near the temperature inversion for models DV2 
(dotted), DV3 (dashed), and DV4 (dash-dotted) at a t — 30000 s, respectively. The initial distribution is shown by the solid line. 



tial model this boundary is located at r = 9.2 10 8 cm and corre- 
sponds to a discontinuous change in the distribution of elements 
(Fig.[3]l, which in stellar evolution models results from the as- 
sumed instantaneous mixing. In such models all species in the 
convectively unstable region are mixed instantaneously across 
the whole convection zone, while the regions which are assumed 
to be convectively stable do not experience any mixing at all. 

The distribution of 12 C at the end of our simulation DV4 



is depicted in Fig. 10 Compared to the initial model there is a 
clear shift of the carbon discontinuity at the outer edge of the 
convection zone to a larger radius (r = 9.7 10 s cm). In hydro- 
dynamic simulations the gas overshoots naturally from the con- 
vectively unstable to the formally convectively stable region be- 
cause of its inertia. At the boundaries of the convection zone the 
overshooting seems always to destroy the stability according to 
the Schwarzschild criterium transforming the originally convec- 
tively stable region into a convectively unstable one. This allows 
the boundary to propagate further when a subsequent load of gas 
will try to overshoot at a later time. We have estimated the prop- 
agation speed of the outer boundary of the convection zone to 
be about ~ 14 m s . With a propagation speed of this order the 
convection zone would reach the hydrogen rich layers surround- 
ing the helium core at a radius r = 1.9 10 9 cm and trigger a hy- 
drogen injection flash (Schlatfl et al.|200l) within just 10 days. 
Expected hydrodynamic phenomena due to the extra hydrogen 
mixing into the helium burning shell via such an extended con- 
vection zone could alter the structure of the star significantly. 
Moreover, additional nucleosynthesis could be triggered, since 
the hydrogen entrainment will result in the production of neu- 
trons and possibly also of some s-process elements. The hydro- 
gen injection flash in Pop I stars is in contradiction to the canon- 
ical scenario since stellar evolutionary models fail to inject hy- 
drogen to the helium core during the core helium flash, unless 
their metallicity is close to zero (Fu jimoto et al.|19 90). 

Since the turbulent entrainment at the inner convective 
boundary involved just three radial grid zones over the longest 
simulations we performed, the estimated propagation velocity 
must be taken with care and be considered as an order of mag- 



nitude estimate. The turbulent entrainment at the outer convec- 
tive boundary involved eighteen numerical zones in radial di- 
mension, therefore the estimated propagation velocity has higher 
confidence level, but still it should be taken as a rough number. 



5.1.3. Two-dimensional models with different resolution 

We find only minor differences between the properties of model 
DV4 and those of the corresponding models computed with a 
different grid resolution. 

First, the initial mapping process leads to different interpo- 
lation errors for different grid resolutions. However, the major 
source of discrepancy in this phase of the calculation is the sta- 
bilization itself. The iterative procedure which minimizes the 
numerical fluxes across zone boundaries (in order to keep the 
model in hydrostatic equilibrium) tends to decrease the temper- 
ature stronger in models with lower resolution. 

Another source of discrepancy is caused by the numerical 
diffusion which is obviously larger in models with lower resolu- 
tion. Therefore, model DV2 suffers more from numerical diffu- 
sion than model DV3 or DV4, which is evident from Figure 1 1 
The temperature inversion, which is almost discontinuous at the 
beginning, gets smoothed out faster in model DV2. Note, that the 
temperature inversion is situated at smaller radii for models with 
higher resolution, since the typical flow velocities are higher in 
better resolved models (Tab.[3 I, i.e., the turbulent entrainment is 
more effective, and the temperature inversion propagates with 
higher speed. 

Nevertheless, models DV3 and DV4 seem to be well re- 
solved since their mutual differences are minor. The temporal 
evolution of their total nuclear energy production rate, for in- 
stance, overlaps almost perfectly (Fig.[TT|i. The temperature fluc- 
tuations in the two-dimensional models are suppressed stronger 
in the better resolved models. Contrary to Dearborn et al. (2006), 
the more intense temperature fluctuations occurring in models 
that we have calculated with grid resolutions even lower than 
that of model DV2, did not lead to an explosion. 
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6. Summary 

We have presented one and two-dimensional {i.e., axisymmetric) 
hydrodynamic simulations of the core helium flash near its peak 
covering about eight hours of evolution time. We find no hydro- 
dynamic events which deviate significantly from the prediction 
of stellar evolutionary calculations. After an initial adjustment 
phase the 2D models reach a quasi-steady state where the tem- 
perature and nuclear energy production rate are only slowly in- 
creasing. 

Convection plays a crucial role in keeping the star in hy- 
drostatic equilibrium. Based on our two-dimensional simulation 
with the highest grid resolution (model DV4), convection fol- 
lows approximately the predictions of mixing length theory, al- 
though the temperature gradient of our dynamically evolved 2D 
models deviates by about 1 % from that of the initial model which 
is obtained from (ID) stellar evolutionary calculations. The max- 
imum temperature (T) max in out best resolved model DV4 rises 
with a rate of about 40 K s _1 , which is about 60% smaller than 
the rate predicted by stellar evolutionary calculations. The mean 
convective velocity exceeds the velocities predicted by mixing 
length theory by up to factor of four. 

During the early 2D dynamic evolution the size of the con- 
vective region does not deviate from that of the initial (hydro- 
static) model. However, after a stable convective pattern is es- 
tablished, our 2D simulations show that the convective flow, con- 
sisting of four quasi-stationary large scale (~ 40 degrees angular 
width) vortices, starts to push the inner and the outer boundary 
of the convection zone as determined by the Schwarzschild sta- 
bility criterium towards the center of the star, and towards the 
stellar surface, respectively. This results in a rapid growth of the 
radial extent of the convection zone on dynamic timescales. 

Our 2D simulations further suggest that it is unlikely that the 
core helium flash is followed by subsequent core helium mini- 
flashes, which are observed in (ID) stellar evolutionary calcula- 
tions, since the inner convective boundary could reach the center 
of the core in less than one month. On the other hand, the injec- 
tion of hydrogen from the stellar envelope into the helium core 
is likely to happen within just 10 days, which is in contradiction 
to the predictions of the canonical evolution of low-mass Pop I 
stars. 

As our 2D axisymmetric simulations probably cannot prop- 
erly capture the intrinsically three-dimensional character of the 
convective flow, we have started to perform also 3D simulations 
of the core helium flash. In addition, we plan to extend our 2D 
simulations to time intervals of several days instead of hours. 
The results of these long-term 2D simulations and of the first 
well resolved 3D simulations of the core helium flash will be 
presented in due time elsewhere. 



1995), one integrates the hydrodynamic equation of energy con- 
servation 



d,(pe) + di(vj(pe + p) - vjLij - Kd{T) = -pv,3,<I> , 

i,j =1,2,3 

(A.l) 

(with e — s + V;V;/2 being the specific total energy density) over 
angular coordinates (6, <p), and separates both the specific en- 
thalpy (e + pip) and the kinetic energy (v,v,/2) into a horizontal 
mean and a perturbation (f = f + /'). This results in 



d,E + d r (F c + F K + F R + F v + F E ) = 



(A.2) 



withQ 
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pe r 2 dQ 

v r p- \s+ -) r 2 dQ 



r 2 dQ , i = 1,2,3 



Kd r T r 2 dQ. 



r z dQ , j = 1,2,3 



Fe — 47IT v r p 



e + — + -VjVi + d r O 
p 2 



(A3) 
(A.4) 
(A.5) 
(A.6) 
(A.7) 

(A.8) 



Here, the various terms F; give the total energy transported per 
unit time across a sphere by different physical processes. They 
are the convective (or enthalpy) flux, F c , the flux of kinetic en- 
ergy, F k, the flux by heat conduction and radiation, Fr, and the 
viscous flux, Fy- Finally, Fe, collects all terms causing a spher- 
ical mass flow, i.e., the model's expansion or contraction, while 
Fc and Fk rest on deviations from this mean energy flow (vor- 
tices). The latter are the major contributors to the heat transport 
by convection, while F v is usually negligibly small. 

In a similar way one can also formulate a conservation equa- 
tion for the mean horizontal kinetic energy, which provides fur- 
ther insight into the effects of convective motions. Using the 
other hydrodynamic equations 
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Appendix A: Energy fluxes 

An analysis of the vertical energy transport allows for conclu- 
sions about the importance of the different physical processes 
occurring in the convection zone. To separate the various con- 
tributions to the total energy flux (Hurlburt et al. 1986] [Achatz | 



d,(p) + diQm) = (A.9) 

d,(pvi) + djidijp + pv iV j - Y.ji) = -p<9,<D , (A. 10) 

i,j= 1,2,3 (A. 11) 

and <9,(pv,v ; /2) = v,d,(pv,) - v,v,<9,p/2, one finds 

d,E K + d r (F K + F P + F V + F e ,k) = P a + Pp + Pv + Pe.k (A. 12) 



1 The gravitational potential O was assumed to be constant for sim- 
plicity. 
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With Fk and Fy as introduced above, one obtains 
E K : 



r 2 dQ 



F P = 



§ 2 ViV ' 
- <T) v r p r 2 dQ 



Fe,k = Anr^-\ P - + ^ 



Pp 



p'djVj r 2 d£l 



(A.13) 

(A. 14) 

(A.15) 

(A. 16) 

(A. 17) 

(A.18) 
(A. 19) 



Pe,k = 4nr 2 -(pdiVi-Vrpd,-®) , i= 1,2,3 

where the f,- are source or sink terms of the kinetic energy. 
They are separated into the effect of buoyancy forces (Pa), fric- 
tion forces (Py), and the work due to density fluctuations (Pp, 
volume changes). By analyzing the various P, one can deter- 
mine what brakes or accelerates convective motions. The acous- 
tic flux, Fp, describes the vertical transport of density fluctua- 
tions. Fe,k and Pe,k describe the effect of expansion (volume 
work, and work against the gravitational potential), similar to 
F E in Eq. (Oj). 
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